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Solvable model of a phase oscillator network on a circle with infinite-range 
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We describe a solvable model of a phase oscillator network on a circle with infinite-range Mexican- 
hat-type interaction. We derive self-consistent equations of the order parameters and obtain three 
non-trivial solutions characterized by the rotation number. We also derive relevant characteristics 
such as the location-dependent distributions of the resultant frequencies of desynchronized oscilla- 
tors. Simulation results closely agree with the theoretical ones. 
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It is ubiquitously observed in nature that a system 
composed of many active elements exhibits collective be- 
havior as a whole. A typical example is the synchro- 
nization of populations of oscillators, e.g., simultaneous 
emission of light by fireflies, the rhythm of the heart com- 
posed of a population of cardiac muscle cells, and circa- 
dian rhythms [i|, Q • 

Pioneering studies on such behavior were done by Win- 
free [3| and Kuramoto [J]. In particular, Kuramoto re- 
garded synchronization as a phase transition and de- 
scribed a prototype model of the phase transition in non- 
equilibrium systems. The model, today called the "Ku- 
ramoto model," is a coupled oscillator system in which 
an oscillator interacts with all other oscillators with the 
same strength. Each oscillator has its own natural fre- 
quency, but its amplitude is constant and the state vari- 
able is its phase. In general, when nonlinear dynamical 
systems with stable limit cycle oscillators are weakly cou- 
pled, the whole system can be described by the phases 
of the oscillators, and the dynamical equation is reduced 
to the evolution equation for phases [4!|. The Kuramoto 
model was used to analytically prove for the first time 
that, as the interaction strength increases from zero, a 
phase transition occurs, from a desynchronized state in 
which each oscillator independently oscillates with its 
own frequency to a synchronized state in which a large 
number of oscillators oscillate with the same frequency 
Q. Since Kuramoto 's analysis of globally coupled os- 
cillators, oscillator networks with short-range and with 
intermediate-range interactions have been studied [6|. 
Oscillators with global and random interactions Q and 
with sparse and random interactions 18] have also been 
studied. Furthermore, the stability of the solutions with 
the Kuramoto model has been studied [y,[i3|- A review of 
the Kuramoto model and its extensions, such as inclusion 
of a noise term, is available elsewhere [ll[. There have 
also been extensive studies on the statistical and dynam- 
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ical properties of the mean-field XY model (HMF XY 
model) of conservative dynamical systems corresponding 
to oscillaor network models of dissipative dynamical sys- 
tems [il,[il. 

Although there have been many studies on oscilla- 
tor networks, no solvable model defined in a finite- 
dimensional space has yet been introduced. It is greatly 
difficult to study systems with short- and intermediate- 
range interactions analytically, so we cannot help re- 
lying on numerical simulation to study such systems. 
To further advance the study of the synchronization- 
desynchronization transition of active elements, it is quite 
desirable to introduce a solvable model that extends the 
Kuramoto model. 

In this paper, we describe a phase oscillator network on 
a circle, and, to make analytical treatment possible, we 
assume infinite-range interaction and that the strength 
and sign of the interaction between two oscillators de- 
pend on the spatial distance between them. We specifi- 
cally adopt the Mexican-hat-type interaction, which was 
introduced to model the creation of feature extraction 
cells in neuroscience and expresses the properties that a 
firing cell excites nearby cells and inhibits distant cells 
|14| . For an XY model on a circle with this interaction, 
it was found that there exists a peculiar solution, the 
pendulum solution, in which the phases of the XY spins 
do not rotate but oscillate as the locations change on the 
circle [l5| . In the phase oscillator network, we show that 
the self-consistent equations (SCEs) of the order param- 
eters for stationary states and the relevant quantities can 
be exactly derived theoretically, and that there exists a 
pendulum solution as in the XY model. 

Now we explain the phase oscillator network. Let (jji 
and di be the phase and location on the circle of the 
j-th oscillator. We regard the i-th oscillator as a two- 
dimensional vector, Xi = {cos (pi, sin (l)i). We assume 
that oscillators are located uniformly on the circle; that 
is, 9i = i27r/N (i = 0, • • • , A^— 1). The evolution equation 



for the i-th phase is 



— (pi = Wi 

at 



— ^ JySin(<^j -(/),)> 



(1) 



which is derived under the rather general situation de- 
scribed above. Here, uji is the natural frequency and is 
drawn from the probability density g(a;), which is as- 
sumed to be one-humped at w = Wq and symmetric with 
respect to wq. In our numerical simulations, we used a 
Gaussian distribution g(w) with a mean of zero and a 
standard deviation a. 

Let us explain the interaction we use in this paper in 
detail. We impose translational symmetry on J^j, i.e., 
Jij takes the form Jy = J{6i — 9j). Furthermore, we 
assume Jij = Jji. Thus, J{0) is an even function of 9. 
Therefore, the Fourier expansion of J{9) is given by 



J(e) = Jo + Ji cos(6l) + J2 cos(26l) -I- 



(2) 



In this study, we treat the case in which only Jq and Ji 
are non-zero, so the interaction is 



'Ji'. 



Jo + Jicos(6'i -9j), 



(3) 



which has the properties expressed by the Mexican-hat- 
type interaction described above. 
The order parameters are defined as 
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Using these order parameters, we rewrite the evolution 
equation ([!]) as 



dt' 



= cjj + Joi?sin(0 — (j)i) 



+Ji[Rc cos 9i sin(0c — (pi) + Rg sin 9i sin(9s — </>i)]. (4) 

Now we derive the SCEs. Without loss of generality, we 
assume ujq = 0. Since we study stationary states, let 
us assume that amplitudes and phases tend to constant 
values as t tends to infinity. We further rewrite eq. (|4]) 
as 



—(jjj = ujj - Aj sin((/)j - aj). 



(5) 



The following relation is derived from a comparison of 
eqs. (HI and ©: 

Aje'"^ = JqRc'^ + JiiRcCOsOje"^" + Rs sm9je'^^]{6) 

Hereafter, we use 9 to identify each oscillator, so Ag is 
expressed as 

Al = {JoRf + Jl{{RcCos9f -t- {Rssm9f 

+2RcRsCos{ec -6s) sin 6* cos 6*} 

+2JoJiR{Rc cos 9c cos 9 + Rs cos Qs sin 9}, (7) 



where 6c = 6c — 6, 6^ = 6s — 6. Defining ipg = (j)g — ae 
transforms the evolution equation into 



-—■Ipg ^UJg - AgSinipg 

at 



(8) 



From this equation, we can develop a theory by following 
Kuramoto's argument. For the synchronized oscillators 
satisfying \u)g\ < Ag, we obtain the entrained phase ipg 
and the number density of the synchronized oscillators 
with the value of phase ip at location 0, ns{9, ip), as 



rg=Sm-\^), 
Ag 



(9) 



nsi9, i)) ^ g{Agsmip)Ag cosip, \ip\ < -, (10) 

where Sin~^x is the principal value and its range is 
[— — , —]■ For the desynchronized oscillators satisfying 

\ujg\ > Ag, we obtain the solution of differential equation 
(|8]) and the number density of the desynchronized oscil- 
lators with the value of phase ip at location 9, nds{9, ip), 
as 



-iAe ^ 


^UJgt + hiujgt), 


(11) 


Cbe -■ 


= wJl-(j^)^ 

V We 


(12) 


nds{9,ip) -■ 


= -/ dx X g{x) ^ .2 ■ 2 r 

T^ J Ae X'' — Ag sm Ip 


(13) 



where ujg is the resultant frequency and h(t) is a peri- 
odic function of t with period 2tt. Note that the en- 
trained phases and the distribution of resultant frequen- 
cies depend on the oscillator locations, in general. From 
eq. (rT3|) . ndsi9,ip + tt) = nds{9,%p) is derived, and 
/o nds{9,ipi)e^'^dip = follows. Thus, only the synchro- 
nized oscillators contribute to the order parameters: 
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(14) 
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dip— I d6i7is(6',?A)sin6ie^'^+^"«,(16) 
27r Jo 



1 /■'" 
where nsiip) — — / d9ns{9,ip). Substituting the ex- 

27r Jp 
pression for ns{9, ip) into these equations, and after some 
algebra, from the real parts of these equations, we obtain 

R = JoR{l) + Ji{Rcf cosOc + Rs9 coses), (17) 

Re = JoRf cos Qc + Ji{Rca + RsCCos{Qc ~ Qs)}, {18) 
Rs = Joi?5coses-H Ji{i?cCCOs(ec-es)+i?s6},(19) 

a = (cos^ 9), h= {s\T? 9),c= (sin 9 cos 9) , 
f = (cose*), g = (sine*), 

^ .7r/2 .27r 

(B) = - di' / d9g{Ag sin iP)cos'^ip B. 

TT Jo Jo 



These are the SCEs for R,Rc, and Rg. Furthermore, we 
derive the following auxiliary equations from the imagi- 
nary parts of eqs. ([T4|) - ([T6l) : 



JoRf sin Qc 
JoRg sin 8 s 



JiRsCsm{Qc 
JiRcCsm{Qc 






0, 
0. 



(20) 
(21) 



From these equations, the phases of the order parameters 
are completely determined. The detailed results will be 
reported elsewhere. There are four solutions of the SCEs, 
and the y are clas sified on the basis of the values of R and 
Ri = ^i?2 + i?2 as 

P: para magnetic solution, {R, Ri) — (0,0), 
U: uniform solution, {R,Ri) = (+,0), 
S: spinning solution, {R,Ri) = (0,+), 
Pn: pendulum solution, (i?, i?i) = (+,+). 

Now, let us consider the physical meanings of these so- 
lutions. To characterize the solutions further, we define 
the rotation number of a solution. The rotation num- 
ber is the number of rotations of synchronized oscillator 
Xg — {cos (jjg , sm(j)g) around the origin in space X as lo- 
cation 9 changes by 27r. In the P solution, all oscillators 
desynchronize, whereas in the other three solutions, an 
extensive number of oscillators synchronize and their di- 
rections are locked. In the U solution, (j)g randomly takes 

a value in the interval [— — ■ + ©, 77 + ©] irrespective of 

the location of the oscillator, so the rotation number is 
0. In the S solution, <j>g linearly depends on 9, and the 
rotation number is 1. In the Pn solution, (jjg fluctuates 
and the rotation number is 0. See Figs. l(a)-(c). 
We display the phase diagram in the scaled parameter 
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FIG. 1. 9 dependencies of entrained phases ^g. Line plots: 
theory; -I-; simulation {N = 10000, cr = 0.2, Jo = 1.2Jo,c). 
Only 1% of entrained phases are depicted, (a) U solution, 
Ji/Jo = 1.9, (b) S solution, Ji/Jo — 2.1, (c) Pn solution, 
Ji/Jo = 2.1. 

space (Jo/cr, Ji/cr) in Fig. 2(a). 

Let us examine the bifurcation structures. As shown 
in the phase diagram in Fig. 2(a), the S and U solutions 
and the S and Pn solutions can coexist. We display the 
Ji/cr dependencies of order parameters R and i?i with 
Jo/cr fixed to 4 in Fig. 2(b). These results show that 
the unstable Pn solution determines the boundary of the 
coexistent regions of the S and U solutions and that of 
the S and Pn solutions. Furthermore, we note that the 
Pn solution continuously bifurcates from the U solution. 
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FIG. 2. (a) Phase diagram in scaled parameter space. Plot 
points represent simulation results; curves represent theoret- 
ical results. Vertical line represents parameters shown in (b). 
(b) Ji/cr dependencies of order parameters. Jo/cr = 4. Solid 
curves represent stable solutions; dashed curves represent un- 
stable solutions, which have superscript U, e.g., S^. 

Taking into account these observations, we derived the 
formulas for the boundaries of bistable regions by using 
the unstable Pn solution and relevant stable solutions. 
In Fig. 2(a), the theoretically obtained boundaries are 
represented by curves. The theoretical results are in good 
agreement with the simulation results. 

Now, let us examine the physical meanings of the phase 
transitions. There are five boundaries in the phase di- 
agram shown in Fig. 2(a). The transition from the 
P to U phase takes place continuously at Jq = Jo,c = 
2/(5(^0)71'), and this is the same transition as in the Ku- 
ramoto model. The transition from the P to S phase 
takes place continuously at Ji = Ji^c = 2Jo,c- In the P 
phase, the rotation number is not defined while it is 1 in 
the S phase. The transition from the U to Pn phase takes 
place continuously at Ji = 2 Jo. In this case, the rotation 
number in both phases is 0. However, as is shown in Figs. 
5(a) and (c), the directions of two synchronized oscilla- 
tors do not correlate in the U phase but correlate weakly 
in the Pn phase because magnitude Ji of the location- 
dependent interaction is larger after the transition than 
before the transition. At the two bistable region bound- 
aries, the stable S and unstable Pn solutions and the 
stable Pn and unstable Pn solutions merge, and the sta- 
ble S and stable Pn solutions disappear. The unstable 
Pn solution differs from the paired solution of the stable 
Pn solution because the phases of the order parameters 
are different. That is, these are a new type of instability 
that does not exist in the Kuramoto model. 

We show theoretical and numerical results for the Jq 




FIG. 3. Jo dependencies of order parameters. Ji = 2.1Jo. 
Curves represent theory; symbols represent simulation results. 
A*' — 20000, a — 0.2. Averages were taken for 20 samples. 
Solid curve and +: R of Pn solution; dashed curve and x: 
-Ri of Pn solution; dashed dotted curve and square: Ri of S 
solution. Vertical lines are error bars. 
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FIG. 4. Distribution G{uj,6) of resultant frequencies for Pn 
solution. Curve represents theory; x represents simulation 
results {N = 100000, cr = 0.2, Ji/Jo = 2.1, Jo = 1.2Jo,c). (a) 
e = 0.05 X 2-K, (b) e = 0.25 X 27r. 



dependencies of the order parameters in Fig. 3, those for 
the location-dependent resultant frequency distribution 
G(w, 9) for different 6 for the Pn solution in Fig. 4, and 
those for the 9 dependencies of the entrained phases (j)g 
for the U, S, and Pn solutions in Figs. l(a)-(c). The 
agreement between the theoretical and numerical results 
is excellent. To investigate the desynchronized oscilla- 
tors, we constructed a Lorenz plot of time series sin(0i(t)) 
for the Pn solution (Fig. 5). The Lorenz plot is defined 




FIG. 5. Lorenz plot for desynchronized oscillator for Pn solu- 
tion, x: theory; +: simulation (A'' — 10000, a = 0.2, Ji/Jo — 
2.1,Jo = l.2Jo,c). 

as the mapping from the difference At; = t;+i — t; to 
Ai/+i, where ti and t;-|_i are the successive times that 
satisfy cos((0i(t/)) = 1 and cos{(j)i{ti+i)) = 1, respec- 
tively. As shown in Fig. 5, the simulation results are 
scattered in the Lorenz plot. This indicates that the 
trajectory of a desynchronized oscillator behaves chaoti- 
cally even though theoretically it is quasi-periodic. How- 



ever, in most of our numerical results, e.g., for the resul- 
tant frequency distribution, the larger the iV, the better 
the agreement between the theoretical and numerical re- 
sults. Our results suggest that the system behaves quasi- 
periodically in the limit of A^ infinity. 

In summary, we have extended the Kuramoto model, 
which is a prototype model of the synchronization- 
desynchronization phase transition in non-equilibrium 
systems, and have proposed a solvable model of a phase 
oscillator network on a circle with infinite-range Mexican- 
hat-type interaction. We derived two auxiliary equations 
by expressing the order parameters by the number den- 
sity of the oscillators. We used them to analytically de- 
termine the phases of the order parameters, derive self- 
consistent equations, and obtain three non-trivial solu- 
tions that are characterized by the order parameters and 
the rotation numbers of the synchronized oscillators X^s. 
We drew phase diagrams by using formulas for the phase 
boundaries derived using the unstable Pn solution, found 
that the unstable Pn solution differs from the paired so- 
lution of the stable synchronized solution, and the transi- 
tion due to pair annihilation of the solution and the rele- 
vant solution is a new type of instability that does not ex- 
ist in the Kuramoto model. We also analytically obtained 
the location-dependent distribution of the resultant fre- 
quencies and entrained phases and validated the theoret- 
ical results by simulation, except for the chaotic behavior 
of the desynchronized oscillators. Our numerical results 
suggest that the system behaves quasi-pcriodically in the 
limit of N infinity. 

In general, when nonlinear dynamical systems that 
have stable limit cycle oscillators are weakly coupled, 
the whole system can be described by phases of oscil- 
lators, and the dynamical equation is reduced to the evo- 
lution equation for phases with general interaction Jij. 
Therefore, by applying the present method to weakly 
coupled dynamical systems on a circle, we should be 
able to obtain new types of solutions and new types 
of synchronization-desynchronization phase transitions. 
Several such studies are now underway. 



[1] D. S. Saunders, An Introduction to Biological Rhythms 
(Blackie, Glasgow and London, 1977). 

[2] A. T. Cloudsley- Thompson, Biological Clocks-Their 
Function in Nature (Weidenfeld and Nicolson, London, 
1980) 

[3] A. T. Winfree, J. Theor. Biol. 16, 15 (1967). 

[4] Y. Kuramoto, Chemical Oscillations, Waves, and Turbu- 
lence ( Springer- Verlag, 1984). 

[5] Y. Kuramoto, in: Proc. Int. Symp. on Mathematical 
Problems in Theoretical Physics, ed. H. Araki (Springer, 
New York, 1975). 

[6] H. Sakaguchi, S. Shinomoto and Y. Kuramoto, Prog. 
Theor. Phys. Lett. 77, 1005 (1987); Y. Kuramoto and 
H. Nakao, Phys. Rev. Lett. 76, 4352 (1996). 

[7] H. Daido, Phys. Rev. Lett. 68, 1073 (1992). 

[8] J. P. L. Hatchett and T. Uezu, Phys. Rev. E 78, 036106 



(2008); J. P. L. Hatchett and T. Uezu, J. Phys. Soc. Jpn., 

78, No. 2, 024001 (2009). 
[9] R. E. MiroUo and S. H. Strogatz, J. Stat. Phys., 60, 245 

(1990); R. E. MiroUo and S. H. Strogatz, J. Nonlinear 

Sci., 17, 309 (2007). 
[10] H. Chiba and I. Nishikawa, Chaos, 21, 043103 (2011). 
[11] J. A. Acebron, L. L. Bonilla, C. J. Perez Vicente and F. 

Ritort, Rev. Mod. Phys. 77, 137 (2005). 
[12] M. Antoni and S. Ruffo, Phys. Rev. E 52, 2361 (1995). 
[13] A. Campa, T. Dauxois and S. Ruffo, Phys. Rep. 480, 57 

(2009) and papers cited therein. 
[14] D. H. Hubel and T. N. Wiesel, J. Physiol. 160, 106 

(1968); thid 195, 215 (1968). 
[15] T. Kimoto, T. Uezu and M. Okada, J. Phys. Soc. Jpn., 

80, No. 7, 074005 (2011). 



